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Abstract 

We present an approximate analytical expression for the escape rate of time-dependent 
driven stochastic processes with an absorbing boundary such as the driven leaky integrate- 
and-fire model for neural spiking. The novel approximation is based on a discrete state 
Markovian modeling of the full long-time dynamics with time-dependent rates. It is valid 
in a wide parameter regime beyond the restraining limits of weak driving (linear response) 
and/or weak noise. The scheme is carefully tested and yields excellent agreement with three 
different numerical methods based on the Langevin equation, the Fokker-Planck equation 
and an integral equation. 
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1 Introduction 



Although the solution of the stationary and unbounded Ornstein-Uhlenbeck pro- 
cess has been found long ago, it is not yet possible to give an analytic exact ex- 
pression that includes time-dependent driving and absorbing boundaries [1,2]. Yet, 
such processes with a linear restoring force an a periodic driving which termi- 
nate at a prescribed threshold are widely used as models for numerous physical 
effects. Examples range from rupturing experiments on molecules [3] where the 
time-dependence is introduced as linear movement of the absorbing boundary up 
to totally different models like the leaky integrate-and-fire (LIF) model for neuronal 
spiking events [4,5,6,7,8]. The latter is the application we primarily think of in this 
paper. The stochastic variable stands for the cell soma's electric potential x(t) that is 
changing due to a great many incoming signals from other neurons. It is thus cus- 
tomary to employ a diffusion approximation for the stochastic dynamics of x(t). 
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The driven abstract LIF model assumes the non-stationary Langevin dynamics (in 
dimensionless coordinates) 

x(t) = -x(t)+Acos(ut + 4>) + V2DZ(t) (1) 

where the process starts at time t — at x(0) = x and fires when it reaches the 
threshold voltage x = a = 1. is white Gaussian noise. Here, a sinusoidal 
stimulus has been chosen for the sake of convenience. The following analysis may 
easily be extended to general periodic stimuli. The dynamics of the process x(t) is 
equivalently described by a Fokker-Planck (FP) equation for the conditional proba- 
bility density function (PDF) p(x, t \ x , 0) in a time-dependent quadratic potential, 

U (x, t) — (x — A cos(tot + <pfj /2, reading 

!, = L( t ),= A („•(,,.),) + i,£p, (2) 

with the absorbing boundary and initial conditions 

p(a, 1 1 xqi 0) = for all t and rr (3) 
p(x, | x ,0) = 5(x - x ). (4) 

After firing the process immediately restarts at the instantaneous minimum of the 
potential. 

The set of eqs. (1-4) defines our starting point for obtaining the firing statistics of 
this driven neuron model. Our main objective is to develop an accurate analytical 
approximation that avoids certain restrictive assumptions of prior attempts. Those, 
in fact, all involve the use of either of the following limiting approximation schemes 
such as the limit of linear response theory (i.e. a weak stimulus A « 1) [7,9] 
or the limit of asymptotically weak noise [10,11,12,13,14]. Our scheme detailed 
below yields novel analytic and tractable expressions beyond the linear response 
and weak noise limit; as will be demonstrated, this novel scheme indeed provides 
analytical formulae that compare very favorably with precise numerical results of 
the full dynamics in eqs. (1, 2-4). The arguments given for the agreement of the 
first-passage time distribution also hold for the residence-time [15] which is not 
further considered here. 



2 Reduction to a discrete model 

The periodicity of the external driving with the period T = 2-k/uj allows one to rep- 
resent the time-dependent solution p(x, t) of the Fokker Planck equation in terms 
of Floquet eigenfunctions and eigenvalues of the Fokker-Planck operator, vi(x, t) 
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and pLi, respectively, [10,16] 



-—Vi(x, t) + L(t)vi(x, t) = HiVi(x, t), (5) 

where the eigenfunctions are periodic in time, integrable in x from — oo to a, and 
fulfill the absorbing boundary condition at x = a 

v l (a,t) = 0. (6) 

The time-dependent PDF can be written as a weighted sum of the Floquet eigen- 
functions 

p(x, t) =J2 c i v i( x , f ) exp(A*»*) (7) 

i 

where the coefficients q are determined by the initial PDF. Note that because of the 
absorbing boundary condition at x = a the total probability is not conserved and 
therefore all Floquet eigenvalues have a non- vanishing negative real part. 

The first main assumption that we impose concerns the value of the potential at 
the boundary: The minimum of the potential must always belong to the "allowed" 
region left of the threshold, and, moreover, the potential difference between thresh- 
old and minimum, denoted by AU(t), must always be larger than at least a few D, 
i.e. AU(t)/D > 4. This assumption implies an exponential time-scale separation 
between the average time t k in which the threshold is reached from the minimum 
of the potential compared to the time r r of the deterministic relaxation towards the 
potential minimum. In the dimensionless units used here r r = 1. For the Floquet 
spectrum this implies the presence of a large gap between the first eigenvalue fii 
which is of the same order as — l/r« and the higher ones which are of the order 
— 1 or smaller. After a short initial time of the order 1, all contributions from higher 
Floquet eigenvalues can be neglected and only the first one survives: 

p(x,t) ~ vi(x,t) exp(yUit) (8) 

In general, the Floquet eigenfunctions and the corresponding eigenvalues are diffi- 
cult to determine. A formal expansion in terms of the instantaneous eigenfunctions 
ipi(x,t) of L(t) fulfilling 

L{t)il> k {x,t) = \ k {t)i) k {x,t) (9) 
is always possible though not always helpful 

Vi(x, t) = J2 d ik (t)i> k (x, t). (10) 

k 

The periodicity of Vi(x, t) and ip k (x, t) implies that the coefficients d ik {t) also are 
periodic functions of time. Expansion (10), together with the Floquet equation (5), 
yields a coupled set of ordinary differential equations for the coefficients d ik (t) [17] 
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Mt) - (A fc (t) - m) d lk (t) = J2 d u(t){ ^- + Mt),Mt) ), (ii) 



where (p k (x,t) denotes the instantaneous eigenfunction of the backward opera- 
tor L + (t) belonging to the eigenvalue X k (t) 

L+{t)ip k {x,t) = \ k {t)ip k {x,t). (12) 

The eigenfunctions ip k (x, t) and (p k (x, t) constitute a bi-orthogonal set of functions 
that always can be normalized such that 

(<p l (t),Mt)) = 8u. (13) 

Here, the scalar product (/, g) is defined as the integral over the real axis up to the 
threshold: 

(f,g)= f dxf(x)g(x) (14) 

J —oo 

With our second assumption we require that the driving frequency oo is small com- 
pared to the relaxation rate in the parabolic potential. Under this condition, the 
matrix elements {d(pk(t)/dt,ipi(t)) that are proportional to the frequency uo are 
also small and may be neglected to lowest order in the equations for the coeffi- 
cients dik(t) [17]. The resulting equations are uncoupled and readily solved to yield 
with the periodic boundary conditions 

dik(t) w S^exp^dt'X^t') - futj, (15) 

where p\ — ^ J T X\[t) dt follows from the periodicity of du(t). Together with 
eqs. (8) and (10) we obtain for the long-time behavior of the PDF 

p(x,t) wexpQ^cft'Ai^)) ipi(x,t). (16) 

Note, that the first Floquet eigenvalue has canceled. The lowest instantaneous eigen- 
functions ip±(x, t) and (fi(x, t) are related by 

il>i(x, t) = (pi(x, t)p {x, t), (17) 

where 

Po(x,t) oc exp(-U(x,t)/Dy (18) 
For the corresponding eigenvalue Xi(t) we find from (9) 

A / t x = S-oo dx M x , t)L(t)v?i(a:, t)p (x, t) 

An explicit expression, valid for high potential differences, can be given after lin- 
earization of U about a 

(pi(x,t) = 1 - exp((a; - a)U'(a,t)/D) (20) 
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which gives for Xi(t) 



X 1 (t) 



AU(t) l-eif(yjAU(t)/P) 
D l-exp(-AU(t)/D) ' 



(21) 



where erf (z) is the error function. 



The waiting-time probability [18] can be expressed as 




(22) 



Therefore, the eigenvalue Xi(t) coincides with the negative of the time-dependent 
escape rate n(t). 

With the expression (21) for the escape rate we can calculate the property of inter- 
est, namely the PDF for the first-passage time (FPT) of the attracting "integrating" 
state that covers the domain — oo < x{t) < a. The FPT-PDF is given by the nega- 
tive rate of change of the waiting time probability, i.e. 



The quantitative validity of these expressions for an extended parameter regime 
will be checked next. 



3 Numerical analysis 

We implemented three different numerical methods to obtain both the FPT-PDF 
and the rate in order to have a reliable basis for comparison with the analytical ex- 
pression (21). The first method performs explicit time-steps of the Langevin equa- 
tion (1). We used an elaborate technique for the time-integration of the fluctuating 
force For points away from the threshold a it is sufficient to take a normal 
distributed random variable for the displacement due to £(t). Quite the contrary in 
the vicinity of the absorbing boundary. Here, the integral of £(t) rather behaves like 
a Wiener process with absorbing boundary, as illustrated in Fig. 1. The appropriate 
transition distribution, is known analytically as the weighted difference between 
two normal distributions [1] 

p(x i+1 ,t + 5t\xi,t) = Ni(x i+ i,5t\xi,0) - N 2 (x i+ i,5t\xi,Q) 



The multiplication on the right-hand side stands for a logical AND that leads to a 
correction step in the algorithm. First, a new position x i+1 is proposed according 




(23) 




(24) 
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to the normal distribution density N±. With the probability P out (xi + i, Xi, St) the 
trajectory has already crossed the boundary during this time-step St from Xi to x i+ i 
and, therefore, is to be ended. The explicit forms of N± and N 2 give 



-Pout^i+l; St) — 




The same formula has been given by [19] but with a different reasoning. 

In order to get the correctly normalized FPT-PDF g(t) we counted the number of 
trajectories hitting the absorbing boundary within the interval [t, t + St). The FPT- 
PDF is then estimated by this number divided by St and by the total number of 
trajectories. The rate is given by 

K(t)=g(t)/P(t), (26) 




Fig. 1. The transition probabilities p(xi + i, t + St | Xj, t) (black line) and N\(x{ + i, t + 5t \ 
X{,t) (dashed line) from Xi for a single time-step, with and without the absorbing boundary, 
respectively. The vertical line indicates the boundary. 
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Fig. 2. First-passage time density g(t) (upper) and rate K,(t) (lower plots) as functions 
of time. Displayed are all three numerical methods we used for testing (solid lines) 
and the approximation based on (21) (dashed lines). The parameters in the left plots, 
AU(t)/D G [5,8] and uj = 0.05, are chosen to yield a very good approximation of the 
rate by (21). The right plots display extreme parameters, AU(t)/D G [3, 8] and uj = 0.5, 
where no good approximation of the rate can be expected. In both cases <fi = 0. 
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where 1 — P(t) is estimated by the number of trajectories that have escaped up to 
time t, divided by the total number of trajectories. 

For the second numerical method we have solved the FP equation (2) using a 
Chebychev collocation method to reduce the problem to a coupled system of or- 
dinary differential equations [13,20]. This gives P(t) as the integral of p(x, t) from 
— oo to a. The FPT-PDF in the figures is then calculated according to eq. (23), and 
the rate again by (26). 

The third method solves Ricciardi's integral equation for the FPT-PDF and is de- 
tailed in [21,22]. For employing his algorithm the process must be transformed into 
a stationary Ornstein-Uhlenbeck process with a moving absorbing boundary 



S(t) = a- A 



1+uj 2 



cos(o>t + (f)) + uj sin(ujt + 0) — e 



(27) 



All three methods provide practically identical results as can be seen in Figs. 2 and 
3. The results for the FPT-PDF and for the rate all collapse into one single line. 
Differences between the numerical methods, e.g. fluctuations in the histogram of 
the Langevin equation method are visible only in the plots of the relative errors 
(Fig. 3 middle and lower rows). 

Figure 2 shows that the FPT-PDF is extremely well approximated by expression 
(21) for the rate n(t). In the left plots we used quite a high barrier with quite slow 
driving compared to the time-scale r r of the process. Good agreement is thus ex- 
pected. In the right plots we show the situation with extreme parameters. The lower 
barrier height A U m j n /D goes down to 3 where a rate-description is unlikely to suf- 
fice. Moreover, the driving is faster, uj = 0.5. The system cannot follow the driving 
instantaneously, and we find a shift in the maximum of the FPT-PDF to later times. 
Under these conditions it is impressing how good the novel approximation still 
works. 



A more delicate measure for the errors of the approximation are the rate K,(t) itself 
and its relative deviation from the three numerically calculated rates. Both can be 
seen in Fig. 3. The upper row of plots shows the approximation error of the rate 
for the same two parameter sets as in Fig. 2. Especially at the maximum the rate 
is over-estimated. This leads to a faster decay of the FPT-PDF which is scarcely 
visible in Fig. 2. Also, the shift of the maxima (indicated by vertical lines) can be 
observed. It is negligibly small for uj = 0.05 but more pronounced for uj = 0.5. 

In the middle row of Fig. 3 a systematic error of the approximation becomes vis- 
ible. The relative error with respect to the numerical results behaves roughly si- 
nus oidally with a phase- shift of 7r/2 relative to the driving and with an additional 
constant offset. For the instantaneous rate expression (19) to be valid it is neces- 
sary that the driving signal is sufficiently slow. If this assumption is violated then 
a rate can still be defined if the barrier is sufficiently high. But in addition to the 
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leading term du(i) in (10) the higher instantaneous eigenfunctions must be taken 
into account. The coupling to the coefficients d lk (t) is induced by the matrix ele- 
ments {(j)k(t),i})i(t)), see eq. (11), containing a time derivative that introduces non- 
adiabatic corrections to the rate and, consequently, to the statistics of the FPT. 

It is quite astonishing, that the huge relative error in the right middle plot of Fig. 3 
leads to such a good result in Fig. 2. The explanation for this is that the FPT- 
PDF (23) uses the time-integrated rate. Therefore, errors are important only where 
the rate is large. A closer look on the plot shows that around the maxima of the rate 
the relative error is comparably small. Because the errors are linear in time around 
the rate's maxima they cancel out when integrated over time in (23). The same is 




Fig. 3. Comparison of the numerical rates and the novel approximation from eq. (21). The 
respective parameters in the left/right plots are the same as in Fig. 2. Upper plots: The 
numerically determined rates are displayed as solid lines: Langevin equation simulations 
(black); Fokker-Planck equation (blue); Ricciardi's integral equation (green). The theoret- 
ical approximation n(t) from eq. (21) is displayed as the red dashed line. Middle plots: 
Relative error of the approximation n(t) with respect to each numerical rate K a (t) (with 
the same color coding as above). Lower plots: Errors of the numerical rates with respect to 
each other. The thin vertical lines indicate the positions of the numerical rates' maxima. 
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Fig. 4. Relative error and relative time-shift of the rate (21). As the basis of comparison we 
used the rate kric obtained by solving Ricciardi's integral equation [21] at times t > r r 
where it has become periodic. Upper plots: The error relative to kric evaluated at the max- 
ima of Kgjc- Shown are data for (AU max /D — AU m i n /D) G (0.1, 1, 2, 3, 5), from top to 
bottom with the colors (red, green, blue, cyan, magenta), and the phase 4> £ (0, 7r/2, 7r) 
with the symbols (x, O, +). In the left panel the driving is slow, uj = 0.05, in the right 
it is fast, uj = 0.5. Note that the relative error is of the same order of magnitude for slow 
and for fast driving. A dependence on the phase cannot be observed. Lower plots: The 
difference of the maxima's position of kric and rate (21) in units of the period T, again for 
uj = 0.05 (left panel) and uj = 0.5 (right panel). Color and symbol codings are the same as 
above. 
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Fig. 5. Relative error and relative time-shift of the maxima as a function of uj. For all 
data points AU m [ n /D = 5. The color and symbol codings for AU max /D and <fi are the 
same as in Fig. 4. For fast driving the relative error starts depending on the interval length 
of AU(t)/D. 
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valid for the residence time whose PDF also contains integrals of the rate [15,23]. 

Figure 4 shows this relative error of n(t) at the maxima of the numerically obtained 
rate as a function of the barrier height. Again, two different driving frequencies are 
given. In both cases the relative error has the same order of magnitude, and thus 
explains why both parameter sets in Fig. 2 yield good approximations. 

Finally, we would like to point the reader's attention to the limitations of the linear 
response approximation. For linear response the parameter ratio A/D needs to be 
small. In our validating example in Fig. 2 (left plots) it takes on the value A/D = 
1.5. Thus, our approximation scheme is valid beyond the linear response limit. 

The time-scale of the driving force is mainly restricted by the relaxation time- 
scale r r and much less by the magnitude of the rate itself. There is no restriction 
on the relative magnitudes of k and T. Instead, both cu and k have to be sufficiently 
small. Fig. 5 indicates that both the relative error and the time-shift of the maxima's 
positions are modest for u < 0.1. 



4 Conclusions 



By reference to a discrete Markovian dynamics for the corresponding full space- 
continuous stochastic process we succeeded in obtaining an analytical approxima- 
tion for the time-dependent escape rate which can be used for calculating first- 
passage time statistics. This result is valid beyond the restraining limits of linear 
response or asymptotically weak noise and of adiabatically slow driving. 

We checked our findings using simulations of the Langevin equation (1) and nu- 
merical solutions of the equivalent FP equation in (2) and of the integral equation 
in [21]. We found an impressive agreement for the first-passage time density and a 
good match for the rate which is the more delicate property for comparison. 

Finally, we note that our method is not restricted to a periodic forcing but applies 
also to arbitrary drive functions. However, in the oscillatory case some of the ap- 
proximation errors cancel out. This leads to useful results in extreme parameter 
regimes where agreement cannot be expected a priori. 

This work has been supported by the Deutsche Forschungsgemeinschaft via project 
HA1517/13-4 and SFB-486, projects A10 and B13. 
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